!python --version
Introduction¶
Dans ce notebook on va essayer de traiter les données des clients d'une banque (ou demandeurs de crédit d'une banque plutôt) avec des outils de Maching Learning, et plus précisément d'apprentissage supervisé.
Apprentissage supervisé¶
C'est un apprentissage où les vraies valeurs de la cible nous sont déjà connues. Alors, nos algorithmes sont comme des enfants, ils apprennent d'une partie des données, on valide leur connaissance sur des observations similairas mais pas vues pendant l'apprentissage et ensuite on teste leur performance sur des observations nouvelles. S'ils ont bien appris, mais pas tout appris par coeur sans comprendre (cf. overfitting) on aura des prédictions assez satisfaisantes pour notre cible.
Les deux grandes classes de l'apprentissage supervisé sont:
Classification - cible discrète, en classes
Regression - cible continue
Cet exercice est un exercice de classification. Nous avons des observations qui sont DEFAULT et d'autres qui ne le sont pas, donc il s'agit d'une classification binaire.
Commençons par importer les modules nécessaires pour les tout premières étapes.
# Modules nécessaires
# Essentiels pour le traitement des données
import numpy as np
import pandas as pd
from math import *
from datetime import datetime as dt
# Affichage
import matplotlib.pyplot as plt
import seaborn as sns
import plotly.express as px
from ipywidgets import widgets
from IPython.display import display
import plotly.graph_objs as go
from plotly.subplots import make_subplots
import plotly.figure_factory as ff
pd.options.plotting.backend = "plotly"
# Fixer la graine
np.random.seed(7)
# Ma palette des couleurs
palette = ['#c95c54', '#00bfdc', '#f7ca28', '#39d4a0', '#9c82c7', '#fa9b4d',
'#c95c54', '#00bfdc', '#f7ca28', '#39d4a0', '#9c82c7', '#fa9b4d']
# Fonction pour affichage en pour cent
def to_percent(df):
return df.apply(lambda x: str(round(x*100,2)) + '%')
Importation des données¶
On lis les données de la base et on affiche des informations importantes sur son contenu.
# Importation des données brutes
df = pd.read_excel("Credit_Cleaned.xlsx")
# Informations sur le dataset
df.info()
# Petit apperçu
df.head()
Nous avons 18 variables explicatives et 1 variable cible.
Variables categorielles:
- Customer_Type
- P_Client
- Educational_Level
- Marital_Status
- Number_Of_Dependent
- Prod_Sub_Category
- Source
- Type_Of_Residence
- Nb_Of_Products
- Prod_Category
- Prod_Closed_Date_NA
Les variables qui contiennent des dates:
- BirthDate
- Customer_Open_Date
- Prod_Decision_Date
- Prod_Closed_Date
Discrètes:
- Years_At_Business
- Years_At_Residence
- Net_Annual_Income
On vérifie s'il y a des données non-renseignées. Puisqu'on a la version nettoyée de ces données, on s'attend à ce qu'il n'y ait pas de données nulles, ce qui est bien le cas:
# Potentielles valeurs Null
df.isnull().sum()
Ensuite, sauvegardons les variables categorielles à part:
categorical_features = ['Customer_Type', 'P_Client', 'Source',
'Nb_Of_Products', 'Prod_Sub_Category',
'Educational_Level', 'Marital_Status',
'Type_Of_Residence', 'Number_Of_Dependant',
'Prod_Category', 'Prod_Closed_Date_NA']
descrete_features = ['BirthDate', 'Customer_Open_Date', 'Years_At_Residence',
'Net_Annual_Income', 'Years_At_Business', 'Prod_Decision_Date',
'Prod_Closed_Date']
print(f"Nombre de variables categorielles: {len(categorical_features)}")
print(f"Nombre de variables discrètes: {len(descrete_features)}")
print(f"Au total: {len(categorical_features) + len(descrete_features)} variables explicatives.")
Ces variables doivent être encodées afin qu'on puisse commencer à tester des différents modèles.
Commençons par encoder la cible afin de pouvoir étudier les liens entre elle et les autres variables.
Dans toutes les observations où la valeur de Y est NO_DEFAULT on mettra 0, et là où il y a DEFAULT on mettra 1.
# Encodage de la cible
print(f"Classes avant encodage:\n{df['Y'].value_counts()}\n")
df['Y'].replace({"NO_DEFAULT" : 0, "DEFAULT" : 1}, inplace = True)
print(f"Classes après encodage:\n{df['Y'].value_counts()}")
# Highly skewed data
to_percent(df['Y'].value_counts() / df['Y'].count())
On remarque que la cible est assez biaisée (skewed en anglais). C'est-à-dire, on a beaucoup plus d'exemples de NO_DEFAULT que de DEFAULT. Plus précisément, environ 92,7% de NO_DEFAULT et 7.3% de DEFAULT.
Donc, si on prédit 0 pour toute observation d'un ensemble test, on aura une justesse (accuracy) très élevée, mais la précision réelle de l'algorithme sera basse, car on ne va jamais obtenir les prédictions importantes, là où il y a un défault, donc les prédictions de la classe 1.
Ceci va jouer un rôle essentiel dans la définition du modèle qu'on va appliquer à ce jeu de données afin d'obtenir un modèle d'une précision satisfaisante.
Analyse et modifications des variables¶
D'abord on crée des fonctions qui nous aideront à mieux analyser les variables explicatives (features) et la cible (target).
Commençons par une fonction qui affiche les histogrammes des variables données en argument. On a la possibilité d'indiquer si les variables qu'on veut afficher sont categorielles ou pas grâce au booléan categorical, et aussi un booléan percent qui sert à indiquer si on veut voir les pourcentages que les groupes représentent de toutes les observations. L'argument nb_cols sert à indiquer le nombre de colonnes de l'affichage de tous ces graphiques, à modifier selon les préférances de l'utilisateur.
def plot_features(features, categorical = True, percent = False, nb_cols = 3):
nb_features = len(features)
histnorm_ = ''
if percent: histnorm = 'percent'
# Make subplots
fig = make_subplots(rows=int(np.ceil(nb_features/nb_cols)), cols=nb_cols,
subplot_titles = features)
for i in range(nb_features):
classes = df[features[i]].value_counts()
classes.sort_values(ascending = False, inplace = True)
nb_classes = classes.shape[0]
if nb_classes < 20:
# Counts if categorical feature
fig.add_trace(go.Histogram(x = df[features[i]],
histnorm = histnorm_,
marker_color = palette,
showlegend = False),
row = int(np.floor(i/nb_cols)) + 1, col = i%nb_cols + 1)
else:
# Normal histogram if descrete/continuous feature
fig.add_trace(go.Histogram(x = df[features[i]],
histnorm = histnorm_,
marker_color = palette[4],
showlegend = False),
row = int(np.floor(i/nb_cols)) + 1, col = i%nb_cols + 1)
fig.update_layout(width = 800, height = 800*nb_features/(nb_cols*3))
fig.show()
Affichons donc, toutes les variables, la cible y incluse, pour avoir un premier apperçu de leurs distributions.
plot_features(df.columns,nb_cols = 3)
Ensuite, on créé une fonction qui nous permet d'afficher l'histogramme d'une variable categorielle à gauche et le pourcentage des observations avec Y=0 et Y=1 selon la catégorie dans laquelle elles se trouvent.
def plot_target_vs_categorical_feature(feature):
classes = df[feature].value_counts()
classes.sort_values(ascending = False, inplace = True)
nb_classes = classes.shape[0]
tmp = df.groupby([feature])['Y'].mean()
fig = make_subplots(rows=1, cols=2,
subplot_titles=(f"Client count per {feature}",
"DEFAULT (Y=1) vs. NO_DEFAULT (Y=0) "))
# Count
fig.add_trace(go.Histogram(x = df[feature],
marker_color = palette,
showlegend = False),
row=1, col=1)
# Percentage of target = 1
fig.add_trace(go.Histogram(x = df[feature],
y = df['Y'],
histfunc = 'avg',
# histnorm = 'percent',
name = "Y = 1",
marker_color = palette[0]),
row = 1, col=2)
fig.add_trace(go.Histogram(x = df[feature],
y = 1 - df['Y'],
histfunc = 'avg',
# histnorm = 'percent',
name = 'Y = 0',
marker_color = palette[1]),
row = 1, col=2)
#fig.update_layout(showlegend = False)
fig.update_layout(barmode = 'stack')
fig.update_layout(width = 800, height = 500)
#fig.update_xaxes(categoryorder = 'total descending')
fig.show()
Affichons les 11 variables categorielles et la répartition des valeurs de la cible en fonction des catégories pour chaque variable categorielle.
for feature in categorical_features:
plot_target_vs_categorical_feature(feature)
On dispose de vriament très peu d'observations de DEFAULT. Il semble que le fait d'être un locataire récent (New_Rent) influe sur le fait que Y soit DEFAULT. La même remarque peut être faite pour la Prod_Category L et pour Prod_Closed_Date_NA = False.
On continue avec les variables discrètes restantes.
# Years_At_Residence
fig = go.Figure()
fig.add_trace(go.Histogram(x = df['Years_At_Residence'],
marker_color = palette[0]))
fig.update_layout(title = "Years_At_Residence")
fig.show()
On doit essayer d'avoir le plus de variables qui ressemblent plutôt à des gaussiennes pour que les modèles (en particulier les modèles linéaires) soient performants. Je décide donc, de modifier cette variable afin d'étaler son histogramme. Je lui applique une transformation avec la fonction racine carrée:
# sqrt(Years_At_Residence)
fig = go.Figure()
fig.add_trace(go.Histogram(x = df['Years_At_Residence'].apply(lambda x: sqrt(x)),
marker_color = palette[6]))
fig.update_layout(title = "Racine carrée de Years_At_Residence")
fig.show()
On voit que les valeurs sont un peu plus centrées et distribuée d'une façon plus normale. Je décide donc de garder cette transformation.
df['Years_At_Residence'] = df['Years_At_Residence'].apply(lambda x: sqrt(x))
# Years_At_Business
fig = go.Figure()
fig.add_trace(go.Histogram(x = df['Years_At_Business'],
marker_color = palette[0]))
fig.update_layout(title = "Years_At_Business")
Cette variable est très déviée également, alors j'applique une transformation racine carrée double. Voici ce que l'on obtient:
# sqrt(sqrt(Years_At_Business))
fig = go.Figure()
fig.add_trace(go.Histogram(x = df['Years_At_Business'].apply(lambda x: sqrt(sqrt(x))),
marker_color = palette[0]))
fig.update_layout(title = "Years_At_Business")
On la garde également.
df['Years_At_Business'] = df['Years_At_Business'].apply(lambda x: sqrt(sqrt(x)))
# Net Annual Income
fig = go.Figure()
fig.add_trace(go.Histogram(x = df['Net_Annual_Income'],
marker_color = palette[0]))
Cette variable est extrêmement déviée, très probablement parce qu'il y a très peu de gens qui ont des gros salaires, comparé au nombre de gens qui sont pauvres, ou payés correctement.
Alors je la modifie en lui appliquant une transformation log.
# log(Net Annual Income)
fig = go.Figure()
fig.add_trace(go.Histogram(x = df['Net_Annual_Income'].apply(lambda x: log(x,10)),
marker_color = palette[0]))
On décide donc, de prendre plutôt le logarithme en base 10 du salaire net annuel pour nos modèles éventuels.
df['Net_Annual_Income'] = df['Net_Annual_Income'].apply(lambda x: log(x,10))
Variables en forme de Dates¶
plot_features(['BirthDate',
'Customer_Open_Date',
'Prod_Decision_Date',
'Prod_Closed_Date'], nb_cols=1)
Créons une nouvelle variable qui va s'appeler duration qui contiendra la période entre Prod_Closed_Date et Customer_Open_Date. Ainsi, on espère se libérer des données très déviées tout en gardant l'information importante qu'elles portent.
duration = df['Prod_Closed_Date'] - df['Customer_Open_Date']
duration = duration.dt.days
duration.name = 'days'
Regardons son histogramme:
fig = go.Figure()
fig.add_trace(go.Histogram(x = duration,
marker_color = palette[0],
showlegend = False))
On remarque qu'elle est assez déviée, alors on applique une transformation log.
fig = go.Figure()
fig.add_trace(go.Histogram(x = duration.apply(lambda x: log(x,10)),
marker_color = palette[0],
showlegend = False))
On décide d'enlever les variables Customer_Open_Date et Prod_Decision_Date et remplacer l'information qu'elles apportent par leur différence, une durée, plus précisément on prend le logarithme en base 10 de ces durées pour que ça ressebmle le plus à une loi gaussienne.
duration = duration.apply(lambda x: log(x,10))
df.drop(columns = ['Customer_Open_Date', 'Prod_Closed_Date'], inplace = True)
df['prod_duration'] = duration
Créons une fonction qui nous permet de convertir les données de format datetime en un format ordinal.
def convert_datetime_featues_to_ordinals(df, features):
for feature in features:
tmp = df[feature].apply(lambda x: x.toordinal())
tmp = (tmp - tmp.mean()) / tmp.std()
df[feature] = tmp.values
convert_datetime_featues_to_ordinals(df, ['BirthDate', 'Prod_Decision_Date'])
df.columns
Voici un petit aperçu de notre data set à présent:
df.head()
Encodage des variables categorielles¶
Pour cela, on va encoder les variables Customer_Type et Prod_Closed_Date à valeurs dans {0,1} et on va utiliser la méthode OneHotEncoding là où c'est nécessaire.
Classes pour lesquelles il faudrait encoder les données:
Variables à 2 classes:
- Y
- Customer_Type
- P_client
- Source
- Prod_Closed_Date
Variables à 3 classes:
- Nb_Of_Products
- Prod_Sub_Category
Variables à 4 ou plus classes:
- Educational_Level (4)
- Marital_Status (5)
- Type_Of_Residence (5)
- Number_Of_Dependant (10)
- Prod_Category (13)
# Nombre de classes pour les variables categorielles
for column in df.columns:
classes = df[column].value_counts()
if classes.count() < 20:
print(f"Variable: {column}\nNombre de classes: {classes.count()}\nClasses:\n{classes}\n\n")
# Encodage de Customer_Type
print(f"Classes avant encodage:\n{df['Customer_Type'].value_counts()}\n")
df['Customer_Type'].replace({"Non Existing Client" : 0, "Existing Client" : 1}, inplace = True)
print(f"Classes après encodage:\n{df['Customer_Type'].value_counts()}")
# Encodage de P_Client
print(f"Classes avant encodage:\n{df['P_Client'].value_counts()}\n")
df['P_Client'].replace({'NP_Client' : 0, 'P_Client' : 1}, inplace = True)
print(f"Classes après encodage:\n{df['P_Client'].value_counts()}")
# Encodage de Prod_Closed_Date_NA
print(f"Classes avant encodage:\n{df['Prod_Closed_Date_NA'].value_counts()}\n")
df["Prod_Closed_Date_NA"].replace({True : 1, False : 0}, inplace = True)
print(f"Classes après encodage:\n{df['Prod_Closed_Date_NA'].value_counts()}")
# Encodage de Source
print(f"Classes avant encodage:\n{df['Source'].value_counts()}\n")
df["Source"].replace({'Sales' : 1, 'Branch' : 0}, inplace = True)
print(f"Classes après encodage:\n{df['Source'].value_counts()}")
Il nous reste à encoder les variables categorielles restantes, à savoir:
- Educational_Level
- Marital_Status
- Prod_Sub_Category
- Type_Of_Residence
- Prod_Category
def OneHotEncoding_(df, features):
for feature in features:
tmp = pd.get_dummies(df[feature])
for column in tmp.columns:
df[column] = tmp[column].values
df.drop(columns = features, inplace = True)
features = ['Educational_Level', 'Marital_Status', 'Prod_Sub_Category',
'Type_Of_Residence', 'Prod_Category']
OneHotEncoding_(df, features)
Une fois nos variables bien choisies, modifiées et encodées, on peut procéder à l'étape suivante et choisir le modèle de prédiction qui conviendrait le mieux.
df.columns
df.dtypes
corr = df.corr()
px.imshow(corr, color_continuous_scale = 'viridis')
On remarque que la plupart des variables ne sont pas fortement correlées entre elles. Les exceptions sont:
- Master/PhD et University corrélées à -0.93
- Married et Single corrélées à -0.93
- E et P corrélées à 0.895
- prod_duration et Customer_Type à 0.65
Je décide de supprimer les variables Master/PhD, Married et P pour cette raison et de continuer à travailler avec le reste des variables.
df.drop(columns = ['P', 'Master/PhD', 'Married'], inplace = True)
Analayse en composantes principales (ACP)¶
Même si l'analyse en composantes principales a été conçue pour les variables continues, afin de réduire la dimensionnalité des variables explicatives, ici on l'utilise simplement d'un le cadre d'une visualisation de l'information contenue dans notre dataset.
C'est-à-dire, on va effectuer une ACP en 2 dimensions, et puis en 3 pour voir si on peut dissocier les clusters NO_DEFAULT{Y=0} et DEFAULT{Y=1} visuellement. Cela pourrait nous donner une idée sur comment ces valeurs sont distribuées dans toute la data set et nous aider dans le choix des variables à choisir pour les modèles.
features = df.columns.drop('Y')
from sklearn.decomposition import PCA
pca = PCA(n_components=2)
components = pca.fit_transform(df[features])
total_var = pca.explained_variance_ratio_.sum()*100
fig = px.scatter(components, x=0, y=1, color=df['Y'],
title=f'Total Explained Variance: {total_var:.2f}%',)
fig.show()
pca = PCA(n_components=3)
components = pca.fit_transform(df[features])
total_var = pca.explained_variance_ratio_.sum() * 100
fig = px.scatter_3d(
components, x=0, y=1, z=2, color=df['Y'],
title=f'Total Explained Variance: {total_var:.2f}%',
labels={'0': 'PC 1', '1': 'PC 2', '2': 'PC 3'}
)
fig.show()
L'ensemble {Y=1} semble être distribué d'une façon quasiment uniforme. Cela veut dire qu'il serait difficile de faire un modèle qui prédira les données avec une bonne précision sur la probabilité de prédire 1 là où la vraie valeur est en effet 1, et c'est cela qu'on priorise dans ce cas-là. Cependant, en prenant toutes les variables dont on dispose, on pourrait arriver à avoir de résultats satisfaisant. Procédons aux modèles.
Modèles d'Apprentissage Supervisé¶
Tout d'abord, séparons notre dataset en features (X) et target (y)
Remarque: Je prends des copies car la dataframe est un objet mutable en python, donc il se comporte comme un passage par référence dans un langage compilé. Puisque je ne veux pas modifier la dataframe initiale, je prends des copies.
X = df.drop(columns = ['Y']).copy()
y = df['Y'].copy()
Ensuite, on va partager notre dataset en trois parties, une sur laquellue on va entraîner (X_train), une qui va nous servir pour la cross-validation (X_test) et une troisième qui va nous servir pour tester les derniers résultats (X_new_data). On va également standardiser les features pour nous assurer que les algorithmes de descente de gradient optimisés (qui se cachent derrière presque tous les modèles de machine learning) puissent converger.
Pour cela, on va se servir d'une utilité du module sklearn:
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.45, random_state = 7)
X_test, X_new_data, y_test, y_new_data = train_test_split(X_test, y_test, test_size=0.25,
random_state = 7)
sc = StandardScaler()
sc.fit(X_train)
X_train = sc.transform(X_train)
X_test = sc.transform(X_test)
X_new_data = sc.transform(X_new_data)
# Pour faire de la cross-validations stratifiée
X_crossval = np.concatenate((X_train, X_test), axis = 0)
y_crossval = np.concatenate((y_train, y_test), axis = 0)
On s'assure que les observations où Y=1 sont assez représentées dans les ensembles qu'on a obtenu:
print(f"Dans X_train: {(y_train.sum() / y_train.shape[0])*100:.3f}% de Y=1")
print(f"Dans X_test: {(y_test.sum() / y_test.shape[0])*100:.3f}% de Y=1")
print(f"Dans X_new_data: {(y_new_data.sum() / y_new_data.shape[0])*100:.3f}% de Y=1")
Ceci correspond bien au pourcentage de ces données dans toute la base de données, donc on peut continuer à tester des modèles d'apprentissage supervisé.
from sklearn import metrics
from sklearn.linear_model import LogisticRegression
model = LogisticRegression(random_state=7)
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
Lorsque la cible est si déviée, il est très important qu'on observe la matrice de confusion, qui est de la forme:
vraies négatifs | faux positifs
faux négatifs | vraies positifs
Les vraies négatifs sont les personnes qui ne présentent pas de défault sur leur crédit et dont on a également prédit qu'il n'y avait pas de défaut.
Les vraies positifs sont les personnes qui présentent de défault sur leur crédit et dont on a également prédit qu'il y avait de défault.
Les faux négatifs sont les personnes pour lesquelles on avait prédit qu'elles ne présenteraient pas de défaut, mais pourtant elles ont le défault.
Les faux_positifs sont les personnes qui ne présentent pas de défault, mais pour lesquelles on avait prédit qu'elles présentaient de défault, à tort.
Regardons ces valeurs pour ce modèle.
true_positives = (((y_test == 1) & (y_pred == 1)) == True).sum()
false_positives = (((y_test == 0) & (y_pred == 1)) == True).sum()
true_negatives = (((y_test == 0) & (y_pred == 0)) == True).sum()
false_negatives = (((y_test == 1) & (y_pred == 0)) == True).sum()
print(f"Vrais positifs: {true_positives}")
print(f"Faux positifs: {false_positives}")
print(f"Vrais négatifs: {true_negatives}")
print(f"Faux négatifs: {false_negatives}")
Il existe une fonction dans le module sklearn, dans la classe metrics, qui nous donne la matrice de confusion directement.
from sklearn.metrics import confusion_matrix
cm = confusion_matrix(y_test, y_pred)
print(cm)
Il est pratique de pouvoir bien visualiser cette matrice de confusion:
def plot_confusion_matrix(cm):
z_text = [[f'True Negatives = {cm[0,0]}',
f'False Positives = {cm[0,1]}'],
[f'False Negatives = {cm[1,0]}',
f'True Positives = {cm[1,1]}']]
fig = ff.create_annotated_heatmap(cm,
x=['0','1'], y=['0','1'],
annotation_text=z_text, colorscale='plasma')
# add custom xaxis title
fig.add_annotation(dict(font=dict(color="black",size=15),
x=0.5,
y=-0.15,
showarrow=False,
text="Predicted value",
xref="paper",
yref="paper"))
# add custom yaxis title
fig.add_annotation(dict(font=dict(color="black",size=15),
x=-0.12,
y=0.5,
showarrow=False,
text="Real value",
textangle=-90,
xref="paper",
yref="paper"))
# adjust margins to make room for yaxis title
fig.update_layout(margin=dict(t=50, l=130))
# add colorbar
fig['data'][0]['showscale'] = True
fig.show()
plot_confusion_matrix(cm)
Pour bien étudier la qualité de ce modèle il est important de calculer la precision et le recall. Ces quantités sont définies ainsi:
PRECISION : De toutes les observations pour lesquelles on a prédit un défault (y_pred = 1), quelle partie a vraiment un défaut (y_test = 1)?
RECALL : Pour quelle partie des observations qui ont un défaut a-t-on réussi à prédire le défaut?
Il est assez intuitif à voir que si on prédit 0 pour toute observation du test on aura une précision à 100%, et une justesse de 92.3% (le pourcentage des observations qui ne présentent pas de défault) mais ce modèle n'a aucun intérêt général.
De même, si on prédit 1 pour toute observation, on aura un recall à 100%. Mais, idem, ce modèle est évidemment inutile.
Donc, notre modèle sera bon s'il a une grande précision et un grand recall.
C'est pour cela qu'on utilise le F-score, également appelé $F_1$. En nottant P la précision et R le recall, le $F_1$ est donné par:
Créons une fonction donc, qui nous calculera ces valeurs pour nos données et nos modèles.
def precision_recall_fscore_accuracy(y_test, y_pred):
true_positives = (((y_test == 1) & (y_pred == 1)) == True).sum()
false_positives = (((y_test == 0) & (y_pred == 1)) == True).sum()
true_negatives = (((y_test == 0) & (y_pred == 0)) == True).sum()
false_negatives = (((y_test == 1) & (y_pred == 0)) == True).sum()
print(f"Vrais positifs: {true_positives}")
print(f"Faux positifs: {false_positives}")
print(f"Vrais négatifs: {true_negatives}")
print(f"Faux négatifs: {false_negatives}")
precision = true_positives / (true_positives + false_positives)
recall = true_positives / (true_positives + false_negatives)
print(f"Precision: {precision:.3f}\nRecall: {recall:.3f}")
F1_score = 2 * (precision * recall) / (precision + recall)
print(f"F1 Score: {F1_score:.3f}")
accuracy = (true_positives + true_negatives) / y_test.shape[0]
print(f"Accuracy: {accuracy:.3f}")
return precision, recall, F1_score, accuracy
precision_recall_fscore_accuracy(y_test, y_pred);
Pour augmenter la probabilité de détecter les défauts on va changer le seuil à partir duquel on décide si on prédit 0 ou 1 dans la fonction sigmoid dans la régression logistique. Pour cela, on créé un vecteur de différents seuils et on compare le $F_1$ score pour chacun. A la fin, on choisit celui qui maximise le $F_1$.
from sklearn.metrics import f1_score
thresholds_ = np.linspace(0,1,30)[1:-1]
F1_score_list = []
for threshold in thresholds_:
y_pred_ = (model.predict_proba(X_test)[:,1] >= threshold).astype(bool)
F1_score_list.append(f1_score(y_test, y_pred_))
tmp = pd.DataFrame(F1_score_list, index = thresholds_, columns = ['F1_score'])
tmp.dropna(inplace = True)
best_threshold_ = tmp.iloc[np.argmax(tmp)]
best_threshold_
fig = go.Figure()
fig.add_trace(go.Scatter(x = thresholds_, y = F1_score_list,
marker_color = palette[0]))
fig.update_layout(title = "F1 Score en fonction de différents seuils")
fig.update_xaxes(title = "Threshold")
fig.update_yaxes(title = "F1-Score")
On obtient que le meilleur seuil est environ 0.31. Appliquons cela à nos prédictions
print("Threshold = 0.5")
precision_recall_fscore_accuracy(y_test, y_pred);
print("\nThreshold = 0.31")
y_modified_pred = (model.predict_proba(X_test)[:,1] >= 0.31).astype(int)
precision_recall_fscore_accuracy(y_test, y_modified_pred);
Ceci a bien augmenté le $F_1$ score, tout en diminuant la "justesse" de la prédiction, mais cette dernière n'est pas considérablement baissé non plus. Donc, si on veut avoir meilleure probabilité de reconnaître le plus de vrais défaults possibles (au prix des faux positifs bien évidemment) alors cela serait la stratégie à adopter.
Ensuite traçons la courbe d'apprentissage, calculée sur 10 validations croisées, en utilisant tous les processeurs (parallélisation), qui montre la valeur moyenne du F score sur le train, et sur le test, en fonction de la taille du train.
# Learning Curve
from sklearn.model_selection import learning_curve
train_sizes, train_scores, test_scores = learning_curve(LogisticRegression(),
X, y,
cv = 10,
scoring = 'f1',
n_jobs=-1,
train_sizes=np.linspace(0.01, 1.0, 50))
train_means = np.mean(train_scores, axis=1)
test_means = np.mean(test_scores, axis=1)
fig = go.Figure()
fig.add_trace(go.Scatter(x = train_sizes, y = train_means,
marker_color = palette[1],
name = 'J_train'))
fig.add_trace(go.Scatter(x = train_sizes, y = test_means,
marker_color = palette[0],
name = 'J_cv'))
fig.update_xaxes(title = 'Train size')
fig.update_yaxes(title = 'F_score')
fig.update_layout(title = "Learning curve")
Il semble qu'on a un grand biais, mais ceci est attendu vu que notre cible est très biaisée.
train_sizes, train_scores, test_scores = learning_curve(LogisticRegression(),
X, y,
cv = 10,
scoring = 'roc_auc',
n_jobs=-1,
train_sizes=np.linspace(0.01, 1.0, 50))
train_means = np.mean(train_scores, axis=1)
test_means = np.mean(test_scores, axis=1)
fig = go.Figure()
fig.add_trace(go.Scatter(x = train_sizes, y = train_means,
marker_color = palette[1],
name = 'Train'))
fig.add_trace(go.Scatter(x = train_sizes, y = test_means,
marker_color = palette[0],
name = 'CV'))
fig.update_xaxes(title = 'Train size')
fig.update_yaxes(title = 'ROC_AUC')
fig.update_layout(title = "Learning curve")
Pour remédier au probllème de grand biais, on pourrait faire les choses suivantes:
- ajouter des nouveles features
- augmenter le paramètre de régularisation C
- ajouter des variables polynômiales à partir des variables existantes
On essayera de voir comment on pourrait augmenter C de façon intelligente. Voici une routine qui teste plusieurs valeurs différentes de C et qui calcule le $F_1$ pour chaque modèle appliqué avec chaque valeur de C. Ensuite, on n'aura qu'à choisir le paramètre C qui maximise le $F_1$
# Choix du meilleur paramètre de régularisation C
from sklearn.metrics import f1_score
C_values = np.linspace(0.05,12,20)
f1_scores_train = []
f1_scores_cv = []
for C in C_values:
model = LogisticRegression(C = C, random_state=0).fit(X_train, y_train)
y_pred = model.predict(X_test)
y_pred_train = model.predict(X_train)
f1_scores_cv.append(f1_score(y_test,y_pred, zero_division=1))
f1_scores_train.append(f1_score(y_train,y_pred_train, zero_division=1))
fig = go.Figure()
fig.add_trace(go.Scatter(x = C_values, y = f1_scores_train,
marker_color = palette[1],
name = 'Train'))
fig.add_trace(go.Scatter(x = C_values, y = f1_scores_cv,
marker_color = palette[0],
name = 'CV'))
fig.update_xaxes(title = 'Regularization parameter C')
fig.update_yaxes(title = 'Scores')
fig.update_layout(title = "Learning curve")
Malheureusement, on remarque que la valeur de C n'a pas de grande influence sur le $F_1$, la valeur par défaut (qui est égale à 1) est donc suffisante et on ne changera pas ce paramètre.
Affichons la courbe ROC qui représente le Taux de Vrais Positifs (TVP) (ce qui n'est rien d'autre qu'un synonime pour recall) en fonction du Taux de Faux Positifs (TFP), calculée à partir de la fonction décision et les vraies labels de la cible.
Soit $TVP$ le Taux de Vrais Positifs, $TFP$ le Taux de Faux Positifs, et donc $V$ ou $F$ signifierait Vrais ou Faux, et $P$ et $N$ positifs et négatifs, respectivement. Nous avons:
$$ TVP = \frac{VP}{VP + FN}$$$$ TFP = \frac{FP}{FP + VN}$$Affichons également la courbe Précision/Recall qui représente les valeurs de la précision en fonction des valeurs du recall (selon les seuils: thresholds). L'intégrale de cette courbe est en fait une très bonne mesure de performance pour les modèles de classification: si toutes les prédictions du modèles sont fausses cette valeur sera égale à 0, si toutes les prédictions sont bonnes elle sera égale à 1. Néanmoins, ici il faut se méfier car la cible est très fortement biaisée. L'AUC est invariante aux changements de seuil de choix (threshold). Donc, il faut se méfier dans les cas où il est préférable de maximiser les vrais positifs, quitte à avoir plus de faux positifs (comme cet exercice où il est préférable de trouver tous les défaults potentiels) ou par exemple lorsqu'on veut minimiser le nombre de faux négatifs, quitte à avoir plus de faux positifs (dans le cas de données médicales et de maladies terminales par exemple).
# ROC and Precision/Recall Curve
from sklearn.metrics import roc_curve, auc, precision_recall_curve, average_precision_score
y_score = model.fit(X_train,y_train).decision_function(X_test)
false_positive_rate, true_positive_rate, _ = roc_curve(y_test, y_score)
roc_auc = auc(false_positive_rate, true_positive_rate)
ap = average_precision_score(y_test, y_score)
precision, recall, thresholds = precision_recall_curve(y_test, y_score)
fig = make_subplots(rows = 1, cols = 2,
subplot_titles = [f"ROC curve (AUC = {roc_auc:.2f})",
f"Precision-Recall curve (AP = {ap:.2f})"])
fig.add_trace(go.Scatter(x = false_positive_rate, y = true_positive_rate,
marker_color = palette[0]),
row = 1, col = 1)
fig.add_trace(go.Scatter(x = recall, y = precision,
marker_color = palette[1]),
row = 1, col = 2)
fig.update_layout(showlegend = False)
Les résultats sont satisfaisants, je choisirais quand même de garder la modification du seuil pour maximiser le nombre de défaults trouvés. Le modèle de régression logistique que je garde serait donc:
logr = LogisticRegression(random_state = 7, C = 1)
Modèle n°2: SVM¶
On va appliquer le modèle SVM pour la classification, donc plus connu sous le nom de SVC
from sklearn.svm import SVC
model = SVC(kernel = 'linear', random_state = 3)
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
precision_recall_fscore_accuracy(y_test, y_pred);
Observons la matrice de confusion:
cm = confusion_matrix(y_test, y_pred)
plot_confusion_matrix(cm)
y_score = model.fit(X_train,y_train).decision_function(X_test)
false_positive_rate, true_positive_rate, _ = roc_curve(y_test.values, y_score)
roc_auc = auc(false_positive_rate, true_positive_rate)
ap = average_precision_score(y_test, y_score)
precision, recall, thresholds = precision_recall_curve(y_test, y_score)
fig = make_subplots(rows = 1, cols = 2,
subplot_titles = [f"ROC curve (AUC = {roc_auc:.2f})",
f"Precision-Recall curve (AP = {ap:.2f})"])
fig.add_trace(go.Scatter(x = false_positive_rate, y = true_positive_rate,
marker_color = palette[0]),
row = 1, col = 1)
fig.add_trace(go.Scatter(x = recall, y = precision,
marker_color = palette[1]),
row = 1, col = 2)
fig.update_layout(showlegend = False)
Ce modèle se comporte mieux que la régression logistique sur nos données de validation. De plus, j'avais testé plusieurs kernels posisbles et le kernel linéaire donnait les meilleurs résultats. Essayons de modifier le paramètre de régularization C.
C_values = np.linspace(0.05,30,20)
from sklearn.metrics import roc_auc_score, f1_score
f1_scores_train = []
f1_scores_cv = []
for C in C_values:
model = SVC(C = C, kernel = 'linear', random_state = 3).fit(X_train, y_train)
y_pred = model.predict(X_test)
y_pred_train = model.predict(X_train)
f1_scores_cv.append(f1_score(y_test,y_pred, zero_division=1))
f1_scores_train.append(f1_score(y_train,y_pred_train, zero_division=1))
fig = go.Figure()
fig.add_trace(go.Scatter(x = C_values, y = f1_scores_train,
marker_color = palette[1],
name = 'Train'))
fig.add_trace(go.Scatter(x = C_values, y = f1_scores_cv,
marker_color = palette[0],
name = 'CV'))
fig.update_xaxes(title = 'Regularization parameter C')
fig.update_yaxes(title = 'F Score')
fig.update_layout(title = "F Score in function of Regularization parametar C")
Il semble que la valeur par défaut (C = 1) est un assez bon choix pour la régularisation. Donc, on le laisse ainsi.
En conclusion, ce modèle donne de très bons résultat sur nos données de test.
svc = SVC(kernel = "linear", C = 1, random_state = 3)
Modèle n° 3: RandomForest¶
Le troisième, et dernier modèle qu'on va essayer c'est le modèle de Forêt Aléatoire.
Les hypermarapètres ici, je les ai cherchés intuitivement car j'ai de l'expérience avec des données similaires, autant en quantité qu'en qualité, et le plus important - en biais. Le paramètre bootstrap est mis à False car comme nos données présentent un grand biais en vue de la cible, tirer des observations par la méthode bootstrap nous donnerait encore plus d'observations de NO_DEFAULT et cela n'aidera pas du tout notre modèle.
Voici ce que l'on obtient:
from sklearn.ensemble import RandomForestClassifier
model = RandomForestClassifier(n_estimators = 150,
max_depth = 12,
bootstrap = False,
min_samples_split = 10,
random_state = 7,
n_jobs=-1)
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
precision_recall_fscore_accuracy(y_test, y_pred);
cm = confusion_matrix(y_test, y_pred)
plot_confusion_matrix(cm)
Ce modèle a de très bonnes pérformances sur ces données. Comparable au SVM qui était très pérformant également. Regardons le rapport de performance et également le AUC:
from sklearn.metrics import classification_report
print(classification_report(y_test, y_pred))
from sklearn.metrics import roc_auc_score
roc_auc_score(y_test, y_pred)
Le AUC est plus petit que celui des autres modèles, mais comme mentionné précédemment, cela n'est pas une métrique complètement bien adapté à notre cas.
On sauvegarde ce modèle et on procède à la dernière étape.
rf = RandomForestClassifier(n_estimators = 150,
max_depth = 12,
bootstrap = False,
min_samples_split = 10,
random_state = 7,
n_jobs=-1)
Tests et comparaison des modèles¶
Voici une comparaison des modèles sur les données test qu'on a séparées avant.
# Comparaison des modèles
print("---------------- Logistic Regression -----------------")
logr.fit(X_train, y_train)
logr_pred = logr.predict(X_new_data)
print(classification_report(y_new_data, logr_pred))
print(f"AUC = {roc_auc_score(y_new_data, logr_pred)}")
plot_confusion_matrix(confusion_matrix(y_new_data, logr_pred))
print("----------------------- SVC ---------------------------")
svc.fit(X_train, y_train)
svc_pred = svc.predict(X_new_data)
print(classification_report(y_new_data, svc_pred))
print(f"AUC = {roc_auc_score(y_new_data, svc_pred)}")
plot_confusion_matrix(confusion_matrix(y_new_data, svc_pred))
print("------------------- Random Forest ---------------------")
rf.fit(X_train, y_train)
rf_pred = rf.predict(X_new_data)
print(classification_report(y_new_data, rf_pred))
print(f"AUC = {roc_auc_score(y_new_data, rf_pred)}")
plot_confusion_matrix(confusion_matrix(y_new_data, rf_pred))
On peut conclure que le modèle de Régression Logistique et le SVM à kernel linéaire sont les plus performants pour ces données, en vue du F-score et de la justesse des prédictions en général.
Afin de privilégier la détéction des défaults, appliquons la baisse du seuil comme précédamment vu:
logr_modified_pred = (logr.predict_proba(X_new_data)[:,1] >= 0.345).astype(int)
precision_recall_fscore_accuracy(y_new_data, logr_modified_pred)
print(f"AUC = {roc_auc_score(y_new_data, logr_modified_pred)}")
plot_confusion_matrix(confusion_matrix(y_new_data, logr_modified_pred))
Cela donne de très bons résultats de prédiction des défauts dans les nouvelles données.
Comparaison par validation croisée stratifiée¶
(K-fold stratified validation)¶
Enfin faison une comparaison par validation croisée stratifiée. On choisi cette méthode pour pouvoir préserver le pourcentage des obervations de chaque classe dans chaque fold de la validation croisée.
Afin de mieux visualiser cette façon de séparer les données en train/test, regardez le graphique suivant:
from sklearn.model_selection import StratifiedKFold
def plot_cv_indices(cv, X, y, group, ax, n_splits, lw=10):
# Generate the training/testing visualizations for each CV split
for ii, (tr, tt) in enumerate(cv.split(X=X, y=y, groups=group)):
# Fill in indices with the training/test groups
indices = np.array([np.nan] * len(X))
indices[tt] = 1
indices[tr] = 0
# Visualize the results
ax.scatter(range(len(indices)), [ii + .5] * len(indices),
c=indices, marker='_', lw=lw, cmap=cmap_cv,
vmin=-.2, vmax=1.2)
# Plot the data classes and groups at the end
ax.scatter(range(len(X)), [ii + 1.5] * len(X),
c=y, marker='_', lw=lw, cmap=cmap_data)
ax.scatter(range(len(X)), [ii + 2.5] * len(X),
c=group, marker='_', lw=lw, cmap=cmap_data)
# Formatting
yticklabels = list(range(n_splits)) + ['class', 'group']
ax.set(yticks=np.arange(n_splits+2) + .5, yticklabels=yticklabels,
xlabel='Sample index', ylabel="CV iteration",
ylim=[n_splits+2.2, -.2], xlim=[0, 100])
ax.set_title('{}'.format(type(cv).__name__), fontsize=15)
return ax
cmap_data = plt.cm.Paired
cmap_cv = plt.cm.coolwarm
n_points = 100
# Sample data
X = np.random.randn(100, 10)
percentiles_classes = [.1, .3, .6]
y = np.hstack([[ii] * int(100 * perc) for ii, perc in enumerate(percentiles_classes)])
# Evenly spaced groups repeated once
groups = np.hstack([[ii]*10 for ii in range(10)])
fig, ax = plt.subplots()
n_splits = 7
cv = StratifiedKFold(n_splits)
plot_cv_indices(cv, X, y, groups, ax, n_splits)
Cela nous permettra de générer des sous-ensembles des données qui contiendront les mêmes distributions de la classe Y, ou du moins la plus proche possbile. Cette façon préserve également les dépendances dans l'ordre du dataset: dans le cas où le paramètre Shuffle = False, toutes les observations dans la classe k dans un des test sets étaient contigües dans y, ou bien séparées dans y par des observations de classes différentes que k.
En utilisant cette manière de séparer notre dataset, on procède en écrivant une fonction qui va nous permettre d'afficher la courbe ROC pour chacun des modèles précédents et de voir ses performances.
from sklearn.metrics import auc, plot_roc_curve
# Run classifier with cross-validation and plot ROC curves
def ROC_with_stratified_cross_validation(classifier, nb_splits):
cv = StratifiedKFold(n_splits = nb_splits)
tprs = []
aucs = []
mean_fpr = np.linspace(0, 1, 100)
fig, ax = plt.subplots(figsize = (10,8))
for i, (train, test) in enumerate(cv.split(X_crossval, y_crossval)):
classifier.fit(X_crossval[train], y_crossval[train])
viz = plot_roc_curve(classifier, X_crossval[test], y_crossval[test],
name='ROC fold {}'.format(i),
alpha=0.3, lw=1, ax=ax)
interp_tpr = np.interp(mean_fpr, viz.fpr, viz.tpr)
interp_tpr[0] = 0.0
tprs.append(interp_tpr)
aucs.append(viz.roc_auc)
ax.plot([0, 1], [0, 1], linestyle='--', lw=2, color='r',
label='Chance', alpha=.8)
mean_tpr = np.mean(tprs, axis=0)
mean_tpr[-1] = 1.0
mean_auc = auc(mean_fpr, mean_tpr)
std_auc = np.std(aucs)
ax.plot(mean_fpr, mean_tpr, color='b',
label=r'Mean ROC (AUC = %0.2f $\pm$ %0.2f)' % (mean_auc, std_auc),
lw=2, alpha=.8)
std_tpr = np.std(tprs, axis=0)
tprs_upper = np.minimum(mean_tpr + std_tpr, 1)
tprs_lower = np.maximum(mean_tpr - std_tpr, 0)
ax.fill_between(mean_fpr, tprs_lower, tprs_upper, color='grey', alpha=.2,
label=r'$\pm$ 1 std. dev.')
ax.set(xlim=[-0.05, 1.05], ylim=[-0.05, 1.05],
title="Receiver operating characteristic example")
ax.legend(loc="lower right")
plt.show()
Validation croisée stratifiée pour la Régression Logistique:¶
ROC_with_stratified_cross_validation(logr, 7)
Validation croisée stratifiée pour le SVC:¶
ROC_with_stratified_cross_validation(svc, 7)
Validation croisée stratifiée pour le Random Forest:¶
ROC_with_stratified_cross_validation(rf, 7)
Conclusion¶
D'après la dernière analyse des modèles par validation croisée stratifiée, on peut constater que la Régression Logistique et le Random Forest montrent des résultats très satisfaisants, comme sur les nouvelles données X_new_data. En combinant les différentes métriques utilisées pour leur évaluation, on peut donner une classification des modèles par performance:
- Régression Logistique (avec ou sans seuil modifié)
- Random Forest
- SVC
Mais il est importent de noter que ces modèles se comportent de façon assez similaire sur notre jeu de données et il n'y a pas de grosses différences dans leurs performances.
Il existe plusieurs modèles différents pour l'apprentissage supervisé, et peu importe celui qu'on choisit, l'important c'est de pouvoir bien intérpréter les données. Même si on n'est pas familier avec toutes les variables explicatives en détails, il faudrait toujours garder un esprit intuitif en permanance, avant de ne pas se perdre dans les nombreuses options et paramètres qui existent et de garder un esprit critique tout au long de l'anaylse, la modélisation et la validation.